home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / SeqPups / appsrc / autoseq.src / CTrace.C < prev    next >
Text File  |  1996-07-05  |  23KB  |  690 lines
  1. //    ============================================================================
  2. //    CTrace.C                                                        80 columns
  3. //    Reece Hart    (reece@ibc.wustl.edu)                                tab=4 spaces
  4. //    Washington University School of Medicine, St. Louis, Missouri
  5. //    This source is hereby released to the public domain.  Bug reports, code
  6. //    contributions, and suggestions are appreciated (to email address above).
  7. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  8. //    Please see CTrace.H for a description of the CTrace object.
  9. //    ========================================|===================================
  10. #include    "CTrace.H"
  11. #include    "CPeakList.H"
  12. #include    "RInclude.H"
  13. #include    "RInlines.H"
  14. #include    <stdlib.h>
  15.  
  16. //template<class T>                            // ATT C++ 3.01 doesn't like this
  17. //vrsn CTrace<T>::version = "94.05.03";        // unfortunately
  18.  
  19. //    ============================================================================
  20. //    CTrace (CTrace constructor)
  21. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  22. //    This method is called upon creation of a new instance of CTrace.  We'll
  23. //    just make sure that the class is properly initialized and allocate a trace
  24. //    if we're given a size.
  25. //    ========================================|===================================
  26. template<class T>
  27. CTrace<T>::CTrace(size_t sz) :
  28.     trace(NULL),
  29.     size(0),
  30.     scale(1),
  31.     mean(0),
  32.     variance(0),
  33.     max(0),
  34.     min(0),
  35.     baseline(0),
  36.     leftCutoff(0),
  37.     rightCutoff(0),
  38.     error_flag(noError)
  39.     {
  40.     error_flag = AllocateTrace(sz);            // okay if sz = 0
  41.     }
  42.  
  43. //    ============================================================================
  44. //    ~CTrace (CTrace destructor)
  45. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  46. //    Called just before deletion of this instance of CTrace.  For now, we'll just
  47. //    ensure that the trace has been deallocated.  In the future, we may wish to
  48. //    implement a dirty bit and check that for writing the trace, etc.
  49. //    ========================================|===================================
  50. template<class T>
  51. CTrace<T>::~CTrace(void)
  52.     {
  53.     delete(trace);                            // free allocated space;
  54.     trace = NULL;                            //   and clear the pointer
  55.     }
  56.  
  57. //    ============================================================================
  58. //    AllocateTrace
  59. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  60. //    Allocates a trace of specified size and sets the size data member.
  61. //    ========================================|===================================
  62. template<class T>
  63. CTError
  64. CTrace<T>::AllocateTrace(size_t sz)
  65.     {
  66.     if (0 != sz)
  67.         {
  68.         if ( (trace = new T[sz]) == NULL )    // allocate space
  69.             return (error_flag=memError);    // ... failed
  70.         size = sz;
  71.         rightCutoff = sz - 1;
  72.         }
  73.     return noError;
  74.     }
  75.  
  76. //    ============================================================================
  77. //    CalculateStats
  78. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  79. //    Computes mean, max, min values for entire trace.
  80. //    NOTES
  81. //    !    It is possible that the average will be incorrectly computed for traces
  82. //        which cause the sum variable to overflow.  This was not a concern for
  83. //        my traces.  Overflow will depend on the size of the trace and the 
  84. //        average value of the points.  One method to overcome this limitation is
  85. //        to keep waypoint averages which do not overflow.
  86. //        ie. Many machines have long double (=double) limits in the approximate
  87. //        range [DBL_MIN = 2.2E-308, DBL_MAX = 1.8E+308].  Assuming an average
  88. //        trace value of 1000, then we could (theoretically) store traces 2E+305
  89. //        points long.  However, assuming an average trace value of 2E+308, we
  90. //        can store only 1 point.  Hopefully your problem will fit somewhere in
  91. //        this range.  ;-)
  92. //    ========================================|===================================
  93. template<class T>
  94. CTError
  95. CTrace<T>::CalculateStats(void)
  96.     {
  97.     register tracepos index;
  98.     long double sum;
  99.  
  100.     if ((size == 0) || (trace == NULL))
  101.         return (error_flag=traceEmpty);
  102.  
  103.     // else...
  104.     max = min = trace[leftCutoff];            // initialize max & min vars
  105.     for (sum=0,index=leftCutoff; index<=rightCutoff; index++)
  106.         {
  107.         if (trace[index] > max)    max = trace[index];
  108.         if (trace[index] < min)    min = trace[index];
  109.                                 sum += trace[index];    // see notes
  110.         }
  111.     mean = (double)sum/(rightCutoff-leftCutoff+1);
  112.  
  113.     for (sum=0,index=leftCutoff; index<rightCutoff; index++)
  114.         sum    += pow(trace[index] - mean,2);
  115.     if (rightCutoff-leftCutoff > 2)
  116.         variance = (double)sum/(rightCutoff-leftCutoff);
  117.     else
  118.         variance = 0;
  119.  
  120.     return noError;
  121.     } 
  122.  
  123. //    ============================================================================
  124. //    Derivative
  125. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  126. //    Calculate the derivative of the trace using the slope of the line between
  127. //    point i-1 and point i+1 as an approximation for the derivative at point i.
  128. //    Returns the derivatives in a CTrace class.
  129. //    
  130. //    +    The trace is returned in the reference parameter.  This should be a
  131. //        (empty) pointer to a CTrace class.
  132. //    +    The derivatives of the first and last points cannot be computed by the
  133. //        method above; therefore the derivative is 2 points shorter than the
  134. //         source trace.
  135. //    !    For the above reason, the derivative at point i is at index i-1 of the
  136. //        derivative trace.
  137. //    !    User is responsible for disposing of *deriv.
  138. //    ========================================|===================================
  139. template<class T>
  140. CTError
  141. CTrace<T>::Derivative(CTrace<double>** deriv)
  142.     {
  143.     double*        dt;                            // pointer to the deriv trace
  144.     tracepos    index;                        // index of source point
  145.  
  146.     if (0 == size)
  147.         return rangeError;
  148.  
  149.     *deriv = new CTrace<double>(size-2);
  150.     if (*deriv == NULL)                        // deriv size = original size
  151.         return (error_flag=memError);
  152.  
  153.     dt = (*deriv)->Trace();
  154.  
  155.     for (index=1;index<=size-2;index++)
  156.         // deriv @ i = slope of line between i-1, i+1
  157.         //   but deriv @ i = dt[i-1]
  158.         dt[index-1] = ((double)trace[index+1]-(double)trace[index-1])/2.0;
  159.  
  160.     (*deriv)->LeftCutoff(leftCutoff-1);
  161.     (*deriv)->RightCutoff(rightCutoff-1);
  162.     (*deriv)->CalculateStats();
  163.  
  164.     return noError;
  165.     }
  166.  
  167. //    ============================================================================
  168. //    PickPeakIndices
  169. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  170. //    Picks local maxima by a concave down method.  Returns in a
  171. //    CSequence<tracepos> all indices, i, for which the derivative at the point
  172. //    i-1 is positive and the derivative at i+1 is negative (note that this
  173. //    implies the second derivative is negative).
  174. //
  175. //    The algorithm is currently as follows:
  176. //    1. for a trace t with indices in [0,t_size-1], we generate a derivative,
  177. //        dt, whose values in index positions [0,t_size-3] correspond to the
  178. //        1st derivatives at t[1..t_size-2].  That is, dt[x] is the derivative at
  179. //        t[x+1].
  180. //    2.    loop for i in [1,t_size-2]
  181. //          if t[i] >= MinPeakHeight then
  182. //            if fabs(dt[i-1]) < ZeroThreshold then
  183. //              the derivative here is 'near' zero... call it a peak
  184. //              append i to peak list
  185. //            else
  186. //              if (i<=t-size-3) and
  187. //                 (dt[i-1] > 0 and dt[i] < 0)
  188. //                we've crossed 0... take index of max(t[i],t[i+1])
  189. //                 if t[i] > t[i+1], add i to peak list
  190. //                 else add i+1 to peak list (t[i+1]>=t[i])
  191. //            
  192. //    If a first derivative is within +/- ZeroThreshold, the peak is added without
  193. //    further consideration.
  194. //
  195. //    !    User is responsible for disposing of *peaks (passed as an argument).
  196. //    *    Remember that the derivative at index i of t is dt[i-1] for the reason
  197. //        mentioned in the Derivative method (above), and that we're using 0-based
  198. //        numbering.
  199. //    !    This function works as described, but several extension to a more
  200. //        general peak picking method are planned.  ie.  concave up & concave down
  201. //        maxima/minima windows, 2nd derivative 0-crossings
  202. //    ============================================================================
  203. template<class T>
  204. CTError
  205. CTrace<T>::PickPeakIndices(
  206.     T            MinPeakHeight,                // consider only peaks >= this value
  207.     double        ZeroThreshold,                // (abs(height) <= this) ==> peak
  208.     CSequence<tracepos>** peaks)            // peak indices returned here
  209.     {
  210.     CTrace<double>*    dt=NULL;
  211.     CTError        CTErr=noError;
  212.     bool        CSErr=FALSE;
  213.     tracepos    index=0;
  214.  
  215.     if (0 == size)
  216.         {
  217.         cerr << "CTrace::PeakPick (size_t): Fatal: Empty trace." << endl;
  218.         return (error_flag=traceEmpty);
  219.         }
  220.  
  221.     CTErr = Derivative(&dt);                // make the derivative and
  222.     if (noError != CTErr)                    // check for errors
  223.         {
  224.         cerr << "CTrace::PeakPick (size_T): Fatal: Derivative failed." << endl;
  225.         delete dt;
  226.         return (error_flag=CTErr);
  227.         }
  228.  
  229.     *peaks = new CSequence<tracepos>;        // make a new CSequence to put
  230.     if (NULL == *peaks)                        // the peaks into
  231.         {
  232.         delete dt;
  233.         return (error_flag=memError);
  234.         }
  235.  
  236.     for (index=leftCutoff;index<=rightCutoff-2;index++)    // index in the trace
  237.         {
  238.         if (trace[index] < MinPeakHeight)    // reject heights < MinPeakHeight
  239.             continue;
  240.  
  241.         // Any derivatives within +/- ZeroThreshold are considered to be
  242.         // first derivative zeros; note that this implicitly includes the
  243.         // case where the derivative exactly equals 0.  This test alone would
  244.         // not discriminate between concave up and concave down (although the
  245.         // MinPeakHeight test would bias toward concave down); for this reason,
  246.         // the second & third tests are required
  247.         if ((fabs((*dt)[index-1]) <= ZeroThreshold)        // deriv is within ZT    
  248.             && ((index>2) && ((*dt)[index-2]>0))        // deriv to left > 0
  249.             && ((*dt)[index]<0) )                        // deriv to right < 0
  250.             CSErr = (*peaks)->Append(index);            // implies concave-down
  251.         else
  252.             // points are outside our ZeroThreshold. Now check to see if
  253.             // derivatives at t[index] and t[index+1] straddle zero; if so,
  254.             // we're at a peak and we'll append index or index+1 to our peak
  255.             // list by comparing t[index] & t[index+1].
  256.             if ( ((*dt)[index-1] > 0) && ((*dt)[index] < 0) )
  257.                 if (trace[index] > trace[index+1])
  258.                     // first point is closer to 0
  259.                     CSErr = (*peaks)->Append(index);
  260.                 else
  261.                     // second point is closer to 0
  262.                     CSErr = (*peaks)->Append(index+1);
  263.  
  264.         if (CSErr)                            // append failed...
  265.             return (error_flag=memError);    //   assume memory error
  266.         }
  267.  
  268.     return CTErr;
  269.     }
  270.  
  271. //    ============================================================================
  272. //    PickPeaks
  273. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  274. //    Calls PeakPicksIndices to get a sequence of tracepos's.  This sequence is
  275. //    converted into a sequence of peak records and the bounds and width of the
  276. //    peak are computed.
  277. //    ========================================|===================================
  278. template<class T>
  279. CTError
  280. CTrace<T>::PickPeaks(
  281.     T        MinPeakHeight,
  282.     double    ZeroThreshold)
  283.     {
  284.     CSequence<tracepos>*    peak_positions;
  285.     CTError                    error;
  286.     uint                    index;
  287.     const double            ONE_HALF = 0.5;
  288.     const tracepos            MAX_WIDTH = 0;
  289.  
  290.     // pick peaks and get just a list of their positions
  291.     error = PickPeakIndices(MinPeakHeight,ZeroThreshold,&peak_positions);
  292.     if (noError != error)
  293.         {
  294.         cerr << "CTrace::PeakPick: Error from size_t PeakPick = "
  295.             << (int) error << endl;
  296.         delete peak_positions;
  297.         return (error_flag = error);
  298.         }
  299.  
  300.     // get the maximum bounds on right and left of each peak
  301.     for (index=0;index<peak_positions->Size();index++)
  302.         {
  303.         PeakRec        pr;
  304.         double        l,r;
  305.         CTError        pwerror;
  306.  
  307.         pr.center = (*peak_positions)[index];
  308.         pr.height = trace[pr.center];
  309.         pwerror = PeakBounds(pr.center,(T)(ONE_HALF*(pr.height+baseline)),\
  310.                                                             l,r,MAX_WIDTH);
  311.         if (noError != pwerror)
  312.             {
  313.             cerr << "CTrace::PeakPick: Fatal error from PeakBounds (" << (int)pwerror
  314.                 << ")." << endl;
  315.             delete peak_positions;
  316.             return (error = pwerror);
  317.             }
  318.         pr.leftBound = l;
  319.         pr.rightBound = r;
  320.  
  321.         peaks.Append(pr);
  322.         }
  323.  
  324.     delete peak_positions;                    // we're done with the peak list
  325.  
  326.     // clip left and right bounds to the midpoint between two peaks in cases
  327.     // where the right bound of peak i overlaps with the left bound of i+1.
  328.     // I know you're thinking that I don't need the peaks.Size() if test, but
  329.     // actually I do... it's unsigned and (unsigned long)(0-1) causes problems.
  330.     if (peaks.Size() > 0)
  331.         for (index=0;index<peaks.Size()-1;index++)
  332.             {
  333.             PeakRec&    pr1 = peaks[index];
  334.             PeakRec&    pr2 = peaks[index+1];
  335.     
  336.             if (pr1.rightBound > pr2.leftBound)
  337.                 pr1.rightBound = pr2.leftBound = (pr1.center+pr2.center)/2;
  338.             }
  339.  
  340.     // compute the peak widths
  341.     for (index=0;index<peaks.Size();index++)
  342.         peaks[index].width = peaks[index].rightBound - peaks[index].leftBound;
  343.  
  344.     return noError;
  345.     }
  346.  
  347. //    ============================================================================
  348. //    PeakBounds
  349. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  350. //    Begins a lateral search within CenterOfSearch +/- MaxWidth (inclusive),
  351. //    looking for either:
  352. //        1. the (real) indices of trace data points which exactly equal
  353. //            Elevation; or,
  354. //        2. the approximation of the positions of the intersection determined
  355. //            by linear interpolation between two points which straddle the
  356. //            desired height.
  357. //    The function returns the peak bounds in the formal parameters
  358. //    LeftIntersection and RightIntersection; the actual return value of the
  359. //    function is an error code.
  360. //    
  361. //    If no bound can be found for an intersection, CenterOfSearch is returned for
  362. //    that value.  This may occur if the peak is near the end of the trace or
  363. //    the MaxWidth is too narrow.
  364. //
  365. //    Specifying 0 for MaxWidth causes PeakBounds to search between CenterOfSearch
  366. //    and the limits of the trace.  It is acceptable to provide a MaxWidth which
  367. //    causes one boundary to be out of range; it will be adjusted to the trace
  368. //    limit automatically.  If both are out of range, an error is returned.
  369. //    An error is also returned if the Elevation is outside [min,max]. (min and
  370. //    max are assumed to be correctly computed already.)
  371. //    ========================================|===================================
  372. template<class T>
  373. CTError    
  374. CTrace<T>::PeakBounds(
  375.     tracepos    CenterOfSearch,
  376.     T            Elevation,
  377.     double&        LeftIntersection,
  378.     double&        RightIntersection,
  379.     tracepos    MaxWidth)
  380.     {
  381.     tracepos    index;
  382.     tracepos    left_bound  = (MaxWidth == 0 ?
  383.                                 0 : CenterOfSearch - MaxWidth);
  384.     tracepos    right_bound = (MaxWidth == 0 ?
  385.                                 size-1 : CenterOfSearch + MaxWidth);
  386.  
  387.     //
  388.     // check ranges
  389.     //
  390.     // return an error if:
  391.     // 1. the user starts a search outside the trace limits, or
  392.     // 2. /both/ bounds are outside trace limits, or
  393.     // 3. Elevation < min
  394.     // 4. Elevation > trace value at CenterOfSearch
  395.     if ( (CenterOfSearch < 0) || (CenterOfSearch > size-1) )
  396.         return rangeError;
  397.  
  398.     if ( (left_bound < 0) && (right_bound > size-1) )
  399.         {
  400.         cerr << "PeakBounds: left or right bound error." << endl;
  401.         return rangeError;
  402.         }
  403.  
  404.     if ( (Elevation < min) || (Elevation > trace[CenterOfSearch]) )
  405.         {
  406.         cerr << "PeakBounds: Elevation out of bounds." << endl;
  407.         return rangeError;
  408.         }
  409.  
  410.     //
  411.     // adjust bounds
  412.     // at least 1 bound is reasonable; adjust the other if necessary
  413.     //
  414.     if (left_bound < 0)            left_bound = 0;
  415.     if (right_bound > size-1)    right_bound = size-1;
  416.  
  417.     //
  418.     // set intersections to CenterOfSearch by default; correct values will
  419.     // be determined or this will be returned to indicate failure to find
  420.     // one or both intersections
  421.     //
  422.     LeftIntersection = RightIntersection = CenterOfSearch;
  423.  
  424.     //
  425.     // begin search in [left_bound,CenterOfSearch]
  426.     //
  427.     for (index=CenterOfSearch; index>=left_bound; index--)
  428.         {
  429.         if (trace[index] == Elevation)
  430.             {
  431.             LeftIntersection = index;
  432.             break;
  433.             }
  434.         // else...
  435.         // if index != left_bound, then check at index-1
  436.         if ( (index != left_bound) &&
  437.              ( (trace[index] > Elevation) && (trace[index-1] < Elevation) ) )
  438.             // we straddle the Elevation
  439.             {
  440.             LeftIntersection =\
  441.               Interpolate(trace[index],index,trace[index-1],index-1,Elevation);
  442.             break;
  443.             }
  444.  
  445.         if (0 == index)
  446.             break;                        // otherwise we'll underflow
  447.         }
  448.  
  449.     //
  450.     // begin search in [CenterOfSearch,right_bound]
  451.     //
  452.     for (index = CenterOfSearch; index<=right_bound; index++)
  453.         {
  454.         if (trace[index] == Elevation)
  455.             {
  456.             RightIntersection = index;
  457.             break;
  458.             }
  459.         // else...
  460.         // if index != right_bound, then check at index+1
  461.         if ( (index != right_bound) &&
  462.              ( (trace[index] > Elevation) && (trace[index+1] < Elevation) ) )
  463.             // we straddle the Elevation
  464.             {
  465.             RightIntersection = \
  466.               Interpolate(trace[index],index,trace[index+1],index+1,Elevation);
  467.             break;
  468.             }
  469.         }
  470.  
  471.     return noError;
  472.     }
  473.  
  474. //    ============================================================================
  475. //    Scale
  476. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  477. //    Scales all trace values by a constant value.
  478. //    ========================================|===================================
  479. template<class T>
  480. void
  481. CTrace<T>::Scale(double relative_scale)
  482.     {
  483.     for (tracepos index=0;index<size;index++)    // scale trace values
  484.         trace[index] = (T)(trace[index] * relative_scale);
  485.  
  486.     scale *= relative_scale;                // new /absolute/ scale
  487.  
  488.     CalculateStats();
  489.     }    
  490.  
  491. //    ============================================================================
  492. //    Translate
  493. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  494. //    Adds constant to all trace values
  495. //    ========================================|===================================
  496. template<class T>
  497. void
  498. CTrace<T>::Translate(T bl)
  499.     {
  500.     if (bl == 0)
  501.         return;
  502.  
  503.     for (tracepos index=0;index<size;index++)    // translate trace values
  504.         trace[index] += bl;
  505.  
  506.     mean += bl;                                // correct mean, max, min
  507.     max += bl;
  508.     min += bl;
  509.     }    
  510.  
  511. //    ============================================================================
  512. //    Smooth
  513. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  514. //    This is an extremely simplistic smoothing routine.  It should be embellished
  515. //    to accommodate the abstract case of arbitrary smoothing windows with
  516. //    arbitrary weighting matrix.  The routine currently smoothes by using a
  517. //    weighted average of the 3 points about a particular index, with the special
  518. //    case of the endpoints handled by throwing out the third points and
  519. //    normalizing.
  520. //
  521. //    !    Computes the new statistics through every iteration.  Although this is
  522. //        a potential bottleneck, it was left in place to ensure consistency of
  523. //        the trace statistics.
  524. //    ========================================|===================================
  525. template<class T>
  526. CTError
  527. CTrace<T>::Smooth(void)
  528.     {
  529.     const double _wt[] = { 0.25, 0.50, 0.25 };
  530.     const double* wt = &_wt[1];                // now reference by wt[-1], wt[ 0],
  531.                                             // and wt[+1].
  532.  
  533.     T*        strace;                            // destination trace
  534.  
  535.     if (0 == size)
  536.         return rangeError;
  537.  
  538.     if ( (strace = new T[size]) == NULL )    // allocate space
  539.         return memError;                    // ... failed
  540.  
  541.     // compute the first and last smoothed points before smoothing the rest of
  542.     // the trace
  543.     // we must divide by the sum of the two weights used to normalize
  544.     strace[0] = (T)((wt[ 0]*trace[ 0]+wt[+1]*trace[+1])/(wt[ 0]+wt[+1]));
  545.     strace[1] = (T)((wt[-1]*trace[-1]+wt[ 0]*trace[ 0])/(wt[-1]+wt[ 0]));
  546.  
  547.     // compute smoothed points for everything in between (that is, [0,size-2])
  548.     for (tracepos index=1;index<=size-2;index++)
  549.         // deriv @ i = slope of line between i-1, i+1
  550.         //   but deriv @ i = dt[i-1]
  551.         strace[index] = (T)(    wt[-1]*trace[index-1] \
  552.                             +    wt[ 0]*trace[index  ] \
  553.                             +    wt[+1]*trace[index+1] );
  554.  
  555.     delete trace;                            // delete original data
  556.     trace = strace;                            // & put the smoothed data in place
  557.  
  558.     CalculateStats();
  559.  
  560.     return noError;
  561.     }
  562.  
  563. //    ============================================================================
  564. //    ReadTrace & ReadAsText
  565. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  566. //    ReadTrace reads size bytes from the file ifp (already opened with rb
  567. //    permissions). Space is allocated for the trace, the size instance variable
  568. //    is adjusted, and the trace is read.  Read trace is fast because it reads a
  569. //    block of data (size of block = # data points * sizeof(data type)).
  570. //
  571. //    ReadAsText reads an endl-delimited list of values from an input file stream.
  572. //    ========================================|===================================
  573. template<class T>
  574. CTError
  575. CTrace<T>::ReadTrace(FILE* ifp, size_t sz)
  576.     {
  577.     CTError        err;
  578.  
  579.     if (trace != NULL)                        // ensure that we don't
  580.         return (error_flag=traceNotEmpty);    // overwrite existing trace
  581.  
  582.     //else...
  583.     if (ifp == NULL)                        // check for valid file pointer
  584.         return (error_flag=badFile);
  585.     
  586.     // else...
  587.     err = AllocateTrace(sz);
  588.     if (noError != err)
  589.         return (error_flag=err);
  590.  
  591.     // else...
  592.     if ( fread(trace, sizeof(T), sz, ifp) != sz )    // fread & error check
  593.         {
  594.         free(trace);                        // fread failed to read size
  595.         size=0;
  596.         return (error_flag=ioError);        // bytes.
  597.         }
  598.  
  599.     CalculateStats();
  600.  
  601.     return noError;
  602.     }
  603.  
  604. template<class T>
  605. CTError
  606. CTrace<T>::ReadAsText(ifstream& is, size_t sz)
  607.     {
  608.     tracepos    index;
  609.     CTError        err;
  610.  
  611.     if (NULL != trace)                        // ensure that we don't
  612.         return (error_flag=traceNotEmpty);    // overwrite existing trace
  613.  
  614.     //else...
  615.     if ((NULL == is) || (is.bad()))            // check for valid file pointer
  616.         return (error_flag=badFile);
  617.     
  618.     // else...
  619.     err = AllocateTrace(sz);
  620.     if (noError != err)
  621.         return (error_flag=err);
  622.  
  623.     // else...
  624.     for (index=0;index<sz;index++)
  625.         {
  626.         is >> trace[index] >> ws;
  627.         if (is.bad())
  628.             {
  629.             DeallocateTrace();
  630.             return (error_flag=ioError);
  631.             }
  632.         }
  633.     CalculateStats();
  634.  
  635.     return noError;
  636.     }
  637.  
  638. //    ============================================================================
  639. //    WriteTrace & WriteAsText
  640. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  641. //    Writes the data pointed to by the trace data member.  These are the write
  642. //    analogs of ReadTrace and ReadAsText; see them for more info.
  643. //    ========================================|===================================
  644. template<class T>
  645. CTError
  646. CTrace<T>::WriteTrace(FILE* ofp)
  647.     {
  648.     if (trace == NULL)                        // check for non-empty trace
  649.         return (error_flag=traceEmpty);
  650.  
  651.     //else...
  652.     if (ofp == NULL)                        // check for valid file pointer
  653.         return (error_flag=badFile);
  654.     
  655.     //else...
  656.     tracepos length = rightCutoff - leftCutoff + 1;
  657.     if ( fwrite(&trace[leftCutoff], sizeof(T), (size_t)length, ofp) != length )
  658.         return (error_flag=ioError);        // fwrite & error ck
  659.  
  660.     //else...
  661.     return noError;
  662.     }
  663.  
  664. template<class T>
  665. CTError
  666. CTrace<T>::WriteAsText(ofstream& os)
  667.     {
  668.     tracepos    index;
  669.  
  670.     if (trace == NULL)                        // check for non-empty trace
  671.         return (error_flag=traceEmpty);
  672.  
  673.     //else...
  674.     //if ((NULL == os) || (os.bad()))            // check for valid file pointer
  675.     // dgg - ^^ that isn't proper -- os isn't a pointer but an object
  676.     if ((os.bad()))            // check for valid file pointer
  677.         return (error_flag=badFile);
  678.  
  679.     //else...
  680.     for (index=leftCutoff;index<=rightCutoff;index++)
  681.         {
  682.         os << trace[index] << endl;
  683.         if (os.bad())
  684.             return (error_flag=ioError);
  685.         }
  686.  
  687.     //else...
  688.     return noError;
  689.     }
  690.